Using the Integrate module¶
In this notebook, we demonstrate how to use the pySCATMECH Integrate module. Its main purpose is to calculate integrated scatter from BRDF models or integrated cross section from local BRDF models over a complex geometry.
In the following, we use some symbols and definitions:
\(f_\mathrm{r}\) : bidirectional reflectance distibution function (BRDF)
\(\mathrm{d}\sigma/\mathrm{d}\Omega\) : differential scattering cross section (DSC)
\(\theta\) : polar angle
\(\phi\) : azimuth angle
\(\Omega\) : solid angle
\(\rho\) : reflectance
\(\sigma\) : cross section
\(\theta_\mathrm{i}\) : incident polar angle
When we perform numerical integration of the BRDF over a solid angle to obtain reflectance, we approximate :raw-latex:`\begin{equation} \rho = \int_\Omega f_\mathrm{r}(\theta_\mathrm{i},\theta,\phi)\;\cos\theta\sin\theta \; \mathrm{d}\theta \; \mathrm{d}\phi \end{equation}` with :raw-latex:`\begin{equation} \rho \approx \sum_j f_\mathrm{r}(\theta_\mathrm{i},\theta_j,\phi_j)\; w_j \end{equation}` where \(\theta_j\) and \(\phi_j\) are sampled polar and azimuth angles, and \(w_j\) are weighting factors.
For local models that define the differential scattering cross section (DSC), we approximate the integrated cross section \(\sigma\) as :raw-latex:`\begin{equation} \sigma = \int_\Omega \frac{\mathrm{d}\sigma(\theta_\mathrm{i},\theta,\phi)}{\mathrm{d}\Omega}\;\sin\theta \; \mathrm{d}\theta \; \mathrm{d}\phi \end{equation}` with :raw-latex:`\begin{equation} \rho \approx \sum_j \frac{\mathrm{d}\sigma(\theta_\mathrm{i},\theta_j,\phi_j)}{\mathrm{d}\Omega}\; \frac{w_j}{\cos\theta_j} \end{equation}`
In this notebook, we will be demonstrating how to integrate the BRDF to obtain the reflectance or the DSC to obtain cross section.
We always start by importing the libraries we need. We are going to do
some graphing, so we need matplotlib.pyplot, and since we are
demonstrating it, we will need pySCATMECH.integrate.
%matplotlib notebook
import matplotlib.pyplot as plt
from pySCATMECH.integrate import *
Defining Geometries¶
We can define various solid angles. Let’s start with a hemisphere:
hemi = Hemisphere()
The Integrator class will be used to perform integration. Its
constructor takes a point spacing and a detector geometry.
The Integrator.PlotSamplingPoints() function plots the integration
points in a manner in which the symbol diameter is proportional to the
weighting factor. There are some optional parameters that we will
discuss later.
There are three different types of samplings available; type = 1 is
the default. For type=, points are spaced evenly in angle space. For
type=2, points are spaced evenly in solid angle. For type=3,
points are evenly spaced in projected solid angle.
Here, we show the three different integration types for the hemispherical detector.
Notice that pySCATMECH provides deg, which is \(\pi/180\). Also,
notice that the sampling weights become smaller where the density of
points becomes greater.
integrator1 = Integrator(2*deg, hemi, type=1)
integrator1.PlotSamplingPoints()
integrator2 = Integrator(2*deg, hemi, type=2)
integrator2.PlotSamplingPoints()
integrator3 = Integrator(2*deg, hemi, type=3)
integrator3.PlotSamplingPoints()
<IPython.core.display.Javascript object>
<IPython.core.display.Javascript object>
<IPython.core.display.Javascript object>
Let’s build some more interesting detection geometries. We start by demonstrating some builtin shapes. One can collect radiation over a circular cone. Here, the cone is centered on \(\theta=45^\circ\) and \(\phi=90^\circ\) and has a half angle of \(20^\circ\).
cone = CircularCone(theta=45*deg, phi=90*deg, alpha=20*deg)
integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
Here, we have an elliptical cone. The half angle in the x direction is 20°, the half angle in the y direction is 45°.
cone = EllipticalCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg))
integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
We can rotate the previous elliptical cone by angle \(\gamma\), in this case 45°.
cone = EllipticalCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg), gamma=45*deg)
integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
We can have a rectangular cone. In this case, it is centered on the surface normal (\(\theta=0,\phi=0\)) and has half angles of 20° and 45°.
cone = RectangularCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg))
integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
And we can rotate that rectangular cone by \(\gamma=45^\circ\):
cone = RectangularCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg), gamma=45*deg)
integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
We can also define arbitrary paths in projected cosine space with a list of \((x,y)\) pairs. The pairs must represent consecutive vertices around the perimeter. Pretty complicated geometries can be generated this way.
cone = ProjectedPolygon([(0.5,0.5),
(0.5,-0.5),
(0,-0.2),
(-0.5,-0.5),
(-0.5,0.5),
(0,0.2)])
integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
We can also combine shapes with logical operations & (and), | (or), and ~ (not). Thus, we can place holes to let radiation enter and exit the hemisphere.
hole1 = CircularCone(theta=70*deg, phi=0, alpha=5*deg)
hole2 = CircularCone(theta=70*deg, phi=180*deg, alpha=5*deg)
detector = Hemisphere() & ~(hole1 | hole2)
integrator = Integrator(1*deg, detector)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
Here is a more complex geometry, created by combining different elements.
circle = CircularCone(theta=0, phi=0, alpha=70*deg)
hole1 = CircularCone(theta=70*deg, phi=0, alpha=5*deg)
hole2 = CircularCone(theta=70*deg, phi=180*deg, alpha=5*deg)
hole3 = CircularCone(theta=0, phi=0, alpha=30*deg)
center = CircularCone(theta=0, phi=0, alpha=20*deg)
detector = circle & ~(hole1 | hole2 | hole3) | center
integrator = Integrator(1*deg, detector)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
Yet another example:
slot = RectangularCone(theta=0, phi=0, alpha=(90*deg, 6*deg))
detector3 = ~slot
integrator = Integrator(1*deg, detector3)
integrator.PlotSamplingPoints()
<IPython.core.display.Javascript object>
Integrate the reflectance¶
So, let’s integrate a BRDF model. We start by creating a model with a
set of parameters. We are going to use Microroughness_BRDF_Model
with a gaussian rough surface. Then, we integrate it over the hemisphere
(with a port provided to let the incident radiation enter and specularly
reflected radiation exit) using Reflectance.
# The model parameters as a dictionary...
parameters = {'lambda' : 0.532,
'substrate' : 4.05+0.05j,
'type' : 0,
'psd' : {None : 'Gaussian_PSD_Function',
'sigma' : 0.05,
'length' : 1
}
}
# Define the model and set its parameters...
model = BRDF_Model("Microroughness_BRDF_Model", **parameters)
# Define the incident angle...
thetaIncident = 6*deg
# The detector is the hemisphere with a hole cut out of it to let in the
# incident light and to let out the specular reflection...
detector = Hemisphere() & ~CircularCone(theta=0*deg, alpha=10*deg)
# Get the integration points and show them...
integrator = Integrator(1*deg, detector)
integrator.PlotSamplingPoints()
# Calculate and print the reflectance ...
ref = integrator.Reflectance(model, thetaIncident)
print("Integrated reflectance = %g" % ref)
<IPython.core.display.Javascript object>
Integrated reflectance = 0.239009
We can see what this looks like as a function of a parameter. In the following, we plot the integrated reflectance as a function of the root-mean-square (RMS) roughness. Notice that we space the points out logarithmically. Also, note how we can define the incident polarization.
# Get the RMS roughnesses in logarithmic spacing...
sigmas = np.exp(np.linspace(np.log(0.0001), np.log(.01), 20))
# Initialize a list of reflectances...
refs = []
# s-polarized incident radiation...
spol = StokesVector(1, 1, 0, 0)
for sigma in sigmas:
# Set the RMS roughness...
model.setParameter("psd.sigma", sigma)
# Calculate the reflectance...
ref = integrator.Reflectance(model, thetaIncident, incpol=spol)
# Add it to the list..
refs.append(ref)
# Plot it...
plt.figure()
plt.loglog(sigmas,refs)
plt.xlabel("$\sigma$")
plt.ylabel("Reflectance")
plt.title("Reflectance versus RMS Roughness")
plt.show()
<IPython.core.display.Javascript object>
Integrated Cross Section¶
Models that inherit Local_BRDF_Model calculate the differential
scattering cross section (DSC) for individual particles or localized
defects. We can use IntegrateCrossSection to obtain the scattering
cross section. For this example, we use the
Bobbert_Vlieger_BRDF_Model, which is an accurate theory for the
scattering by a sphere above a surface.
BVparameters = {'lambda' : 0.532,
'substrate' : "silicon",
'type' : 0,
'sphere' : 1.59,
'radius' : 0.05,
'spherecoat' : "No_StackModel",
'stack' : {None : 'SingleFilm_StackModel',
'material' : '(1.5,0)',
'thickness' : '0.0016'},
'delta' : 0} # use defaults for operational parameters
BVmodel = Local_BRDF_Model("Bobbert_Vlieger_BRDF_Model",**BVparameters)
thetaIncident = 70*deg
circle = CircularCone(theta=0, phi=0, alpha=70*deg)
hole1 = CircularCone(theta=70*deg, phi=0, alpha=5*deg)
hole2 = CircularCone(theta=70*deg, phi=180*deg, alpha=5*deg)
hole3 = CircularCone(theta=0*deg, phi=0, alpha=30*deg)
detector = circle & ~(hole1 | hole2 | hole3)
integrator = Integrator(2*deg, detector)
integrator.PlotSamplingPoints()
cs = integrator.CrossSection(BVmodel,thetaIncident,incpol=[1,-1,0,0])
print("Cross section = %g µm^2" % cs)
<IPython.core.display.Javascript object>
Cross section = 0.000183612 µm^2
We can see what the cross section is as a function of sphere diameter.
diameters = np.exp(np.linspace(np.log(0.05),np.log(2),20)) #logarithmic spacing
refp = []
refs = []
spol = Polarization('s')
ppol = Polarization('p')
for d in diameters:
BVmodel.setParameters(radius=d/2)
refs.append( integrator.CrossSection(BVmodel, thetaIncident, incpol=spol) )
refp.append( integrator.CrossSection(BVmodel, thetaIncident, incpol=ppol) )
plt.figure()
plt.loglog(diameters,refs,label="s")
plt.loglog(diameters,refp,label="p")
plt.xlabel("Diameter / $\mathrm{\mu m}$")
plt.ylabel("Cross section / $\mathrm{\mu m}^2$")
plt.title("Cross Section versus Diameter")
plt.legend()
plt.show()
<IPython.core.display.Javascript object>
Visualizing Intensity¶
We can visualize the scatter distribution by using the
PlotSamplingPoints with some arguments. The function DSC returns
an array of differential scattering cross sections for each of the
sampling points and we can pass that to the color argument of
PlotSamplingPoints. Because the default symbol size for displaying
the sample points leaves gaps between the symbols, we can expand their
size to fill those gaps using the expand argument.
hemi = Hemisphere()
integrator = Integrator(1*deg, hemi)
dsc = integrator.DSC(BVmodel, thetaIncident, incpol=ppol)
integrator.PlotSamplingPoints(color=np.log(dsc), expand=3)
<IPython.core.display.Javascript object>
Polarization¶
This package is designed to simulate integrated scatter for directional
incident irradiation. The incident polarization is therefore pretty
straightward to define and has been illustrated above using Stokes
vectors. The detector system, however, may also be polarization
sensitive. We use the sensitivity argument of the SolidAngle to
define that sensitivity. The Sensitivity function differs from the
Polarization function, because a factor of two is usually needed
when a detector is sensitive to a single polarization and not the other.
detector = CircularCone(theta=45*deg, phi=90*deg, alpha=30*deg, sensitivity=Sensitivity('s'))
integrator = Integrator(1*deg, detector)
diameters = np.exp(np.linspace(np.log(0.05), np.log(2),20)) #logarithmic spacing
ref = []
for d in diameters:
BVmodel.setParameters(radius=d/2)
ref.append( integrator.CrossSection(BVmodel, thetaIncident, incpol=Polarization('p')) )
plt.figure()
plt.loglog(diameters, ref)
plt.xlabel("Diameter / $\mathrm{\mu m}$")
plt.ylabel("Cross section / $\mathrm{\mu m}^2$")
plt.title("Cross Section versus Diameter")
plt.show()
<IPython.core.display.Javascript object>
